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ABSTRACT 

Context. Understanding the gas abundance distribution is essential when tracing star formation using molecular line observations. 
Changing density and temperature conditions cause gas to freeze-out onto dust grains, and this needs to be taken into account when 
modeling a collapsing molecular cloud. 

Aims. This study aims to provide a realistic estimate of the CO abundance distribution throughout the collapse of a molecular cloud. 
We provide abundance profiles and synthetic spectral lines which can be compared to observations. 

Methods. We use a 2D hydrodynamical simulation of a collapsing cloud and subsequent formation of a protoplanetary disk as input for 
the chemical calculations. From the resulting abundances, synthetic spectra are calculated using a molecular excitation and radiation 
transfer code. 

Results. We compare three different methods to calculate the abundance of CO. Our models also consider cosmic ray desorption and 
the effects of an increased CO binding energy. The resulting abundance profiles are compared to observations from the literature and 
are found to agree well. 

Conclusions. The resulting abundance profiles agree well with analytic approximations, and the corresponding line fluxes match 
observational data. Our developed method to calculate abundances in hydrodynamical simulations should greatly aid in comparing 
these to observations, and can easily be generalized to include gas-phase reaction networks. 
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1. Introduction 

The study of low-mass Young Stellar Objects (YSOs), from 
early protostellar cores to T Tauri stars surrounded by disks, 
involves observations of either the thermal emission from dust 
grains or molecular line emission. While dust e mission yields 
inform ation about the temperature profile (e.g., iDraine & Led 
Il984l) . and, through measurements of the shape of the Spectral 
Energy Distribution, the evolutionary stage dLada & Wilkin g| 
1984; Adams et al. .1 1 1987f) . spectral lines are the only way to con- 
strain the kinematical properties of YSOs. The interpretation of 
molecular lines are crucial for the understanding of protoplan- 
etary disk formation and the distributi on of angu lar momentum 
during these stages of star formation dEvansI 19991) . However, de- 
riving the gas motion from spectral line profiles is complicated 
by degeneracies between the velocity field topology, inclination 
of the angular momentum axis, optical depth, and geometrical 
effects (iBrinch et al.ll2007h . 

In addition to these effects, molecular abundances may vary 
greatly. The gas abundance is determined by the ongoing chem- 
istry driven by the co-evolving density and temperature distribu- 
tions, and to some extent also the velocity field. The variation in 
molecular abundances throughout a YSO may have a huge im- 
pact on line profiles and therefore it is important to understand 
the chemistry in order to interpret these line profiles. The effect 
of the chemistry is enhanced when observations at higher reso- 
lution are done. If the observations are done at a relatively low 
resolution (>10"), the emission is smeared out over a large por- 
tion of the object and therefore an average constant abundance 
is mostly adequate. However, when going to a higher resolution, 
using sub-millimeter interferometry the emission is aver- 



aged over a smaller part of the object and thus variations in the 
abundance need to be taken into account. 

The problem of determining the abu ndance dis t ributio n has 
been addressed before in the literature. iLee et aL I (120041) took 
a semi-analytical approach using a self-similar collapse of an 
isothermal sphere in which they followed "fluid-elements". In 
that paper, the chemistry was followed in a series of Bonnor- 
Ebert spheres during the pre-stellar phase and in a self-similar 
inside-out collapse during the protostellar phase. This model 
provides a good analytical description of the early stages of star 
formation, but, l acking any rot a tional velocity, no d isk is formed 
in this scenario. lAikawa et all d200lL 120031, [2005) followed the 
chemical evolution of contrac t ing pr e-stellar clouds using an 
approach similar to ILee et all d2004l) but using an isothermal 
cloud collapse only. Their calculations consider only the very 
early stages of the star formation process. |j0rgensen et al.l (12002, 
2005) introduced the "drop model", where chemical depletion 
occurs as a step function when certain temperature and density 
conditions of any underlying model are met. While this approach 
is simple yet quite successful, there may be cases where the 
abundance changes only very gradually so that a step function 
is no l onger a good approximation. More recently, Tsami s~et al.l 
(2008) calculated HCO + , CS, and N 2 H + abundances and spectra 
from a simulation of an inside-out collapsing core. 

In this paper we investigate detailed evolution of the gas- 
phase and solid state abundance of CO in a collapsing, rotating 
core. We base our study on a hydrodynamical simulation, where 
we follow the CO freeze-out and evaporation for a number of 
test particles that flow with the gas. The lay-out of the paper is 
as follows: In Section [2] we present our numerical simulations 



and the equations used to solve the CO freeze-out chemistry. In 
Section[3]we show three different chemical models and how our 
results compare to observations as well as the results of previ- 
ous papers. Sections [4] and [5] hold a discussion and a summary 
respectively. 



2. Tracing chemistry during disk formation 

2. 1 . The physical model 

Instead of using an analytic description of the gravitational 
collapse and subsequent disk formation, we use a grid-based 
2D axis-symmetric hydrodynamical scheme. The code which is 
used is described in lYorke & Bodenheimerl(ll999l) and in partic- 
ular we adopt the J-type model described in that paper. 

The initial conditions for this model are a 1 M isother- 
mal sphere with a temperature of 10 K and a power-law den- 
sity slope of p oc r~ 2 . The cloud has an outer radius of lxlO 15 
m (6667 AU), and an initial solid-body rotational perturbation 
of lxlCT 13 s _1 is given. The model evolves under the action 
of gravity while the temperature is solved for self-consistently 
using an approximate radiation transfer method. Angular mo- 
mentum is t ransferred through artifici al viscosity using an a- 
prescription (IShakura & Sunyaevll973l) . 

The age of the simulated system is described in terms of the 
initial free-fall timescale, 



tff = 



7T 2 R 3 

8GM' 



(1) 



where M is the mass of the cloud and R is the radius. For the 
initial conditions in our particular simulation one free-fall time 
is 1 x 10 5 yr. We let the simulation run for about 2.6 tff at which 
point the simulation terminates due to extreme velocities near 
the center. Figure [T| shows four time snapshots at characteristic 
ages. We use the same four snapshots throughout this paper, at 
times of 0.0, 0.5, 1.5, and 2.5 tff (with the exception of Fig.Q] 
where we use a few snapshots later than t = 0.0 in order for the 
temperature to have evolved slightly from the isothermal initial 
condition). 

The luminosity is given by the sum of the intrinsic stellar 
luminosity and the accretion luminosity 



3 GM*M 

4 ~R~, ' 



(2) 



where M„ is the mass of the star, R* the radius of the star, and M 
the mass flux onto the star. The total luminosity is shown, plotted 
as a function of free-fall time, in Fig. [2] After the initial increase 
in the luminosity, it drops for a while, after which it steeply in- 
creases again and then slowly fal ls of linearly for the remainde r 
of the simulation. As described in lYorke & Bodenheimerl (fl999). 
the smooth increase and the drop around 0.5 tff is dominated by 
the accretion luminosity where low angular momentum material 
is able to fall directly onto the star. The sharp increase in lumi- 
nosity marks the formation of an equilibrium disk and from then 
on, the intrinsic stellar luminosity dominates the total luminos- 
ity. 

The hydrodynamical simulations does not include any chem- 
istry and because of that, we cannot let the state of the molecules, 
whether they are in the gas phase or locked in an ice matrix, 
couple back to the hydrodynamics. We therefore determine the 
abundance by post-processing the result of the hydrodynamical 
calculations. We follow the chemistry by populating the compu- 
tational domain of the first snapshot with trace particles. These 
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Fig. 2. The time evolution of the luminosity of the central source 
in the hydrodynamical simulatio n. This figure compares to Fig. 7 
in lYorke & Bodenheimerid 19991) . 



particles are massless and do not interact with each other in any 
way. The trace particles are positioned evenly in the radial di- 
rection and 9 x 10 5 particles are used. Every trace particle rep- 
resents the local molecular environment and carries informa- 
tion on the state of the molecules (ice or gas). The particles 
are then allowed to follow the flow, predetermined by the hy- 
drodynamical calculations, and the state of the particles is up- 
dated as temperature and density conditions change throughout 
the hydrodynamical simulation. Finally, the state of the trace 
particles are mapped back onto the hydro-grid at each output 
time step, so that a complete history from t=0.0 tff to t=2.5 
tff of the CO abundance as function of R and z is linked to 
the density and temperature. Taking the CO distribution into ac- 
count, we can use each of these time snapshots as input models 
for our 2D line excitation and rad iative transfer code RATRAN 
dHogerheiide & van der T ak 2000) to predict line profiles of CO 
and its optically thin isotopologue C 18 0. 

Before proceeding with the freeze-out calculations, it is in- 
teresting to consider the dynamics of the trace particles in order 
to develop an understanding of the environments that the parti- 
cles traverse. By a simple radial color coding of the particles, 
we can effectively visualize the flow within the hydrodynamical 
simulation and track how material is accreted onto the disk. Four 
snapshots of the particle distribution are shown in Fig. [3] where 
the shading of the particles refers to their initial distance from 
the center. The disk is seen to be layered vertically rather than 
radially which is somewhat counterintuitive. The vertical layer- 
ing arises because material joins the disk from above. As the 
disk spreads viscously due to conservation of angular momen- 
tum, a shock front propagates outwards, pushing aside infalling 
material. The only possible path of the particles is therefore up- 
wards and over the disk lobe until they can rain down onto the 
disk at smaller radii. This rolling motion, is also reflected in the 
strong vertical mixing in the outer disk. However, this behavior 
may be an artifact related to the 2D nature of the hydrodynam- 
ical code. In a 3D simulation, particles would also be able to 
dissipate in the azimuthal direction, which would likely change 
their motions. 



2.2. Freeze-out and evaporation 

In general, molecular gas abundances depend on reaction rates 
and phase transitions. However, under the physical conditions 
present in our model, only the latter has any significant impact 
on the CO abundance. The state of CO is governed by rate equa- 
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Fig. 1. The density (left panel) and temperature (right panel) contours of four time snapshots of the hydrodynamical simulation. 
Time progresses from left to right, top to bottom at t = (0.0(6), 0.5, 1.5, 2.5} iff. In the right panel, contour lines refer to the dust 
temperature, while the color scale denote the gas temperature. 
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Fig. 3. Trace particle motion in the hydrodynamical simulation. 
The particles are color coded according to their initial radial po- 
sition. Contour lines describe the density. The vertical color bar 
shows the particles initial distance from the center. The four pan- 
els (from the top-left to bottom-right) show the particle positions 
att=0.0,0.5, 1.5, and 2.5 t//. 



tions. The two main mechanisms are depletion, where gas phase 
molecules freeze-out onto grains, and desorption, where solid 
state molecules evaporate into the gas phase. 



The d epletion rate A used i n this paper is given by the equa- 
tion (e.g.,|ChirnieiitaD|200l, 



2 8kT g 
A — na - 



\ miico 



(3) 



where T g is the gas temperature, mco is the mass of a CO 
molecule, and n gra j„ is the grain number density, assuming a 
grain abundance of 1 .33x 10~ 12 relative to the hydrogen nucleon 
density. For the mean grain size a we adopt the value of 0. 1pm. 
We also assume a sticking probability of unity. 

A similar equation can be written for the desorp- 
tio n of molecules. The th ermal desorption rate is given 
bv lWatson & Salpeteild 1972b . 
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where ky> is Boltzmann's constant, T c \ is the dust tempera- 
ture, vyQ_js_Jhe_vibrational frequency of X in its binding 
site dHasegawa et al.ll 1992b . and E b (X) is the binding energy of 
species X onto the dust grain. In this paper we use a default bind- 
ing energy of 960 K, but other values, corresponding to mantles 
which are not composed of pure CO ice, are explored as well in 
Sect. El 

A second desorption mechanism which is taken into account 
is cosmic r ay induced desorption. Th is mechanism was first 
proposed by IWatson & Salpeterl dl972b . Energetic nuclei might 
eject molecules from grain surfaces by either raising the tem- 
perature of the entire gr ain or by spot heating near the impact 
site. The formulation of lHasegawa & Herbsa dl993b is used to 
calculate the cosmic ray desorption rate, 



= 3.16 x lO" 19 ^ 



T d =10K 



(5) 



Cosmic ray desorption can be added to the thermal desorption, 
effectively preventing a depletion of 100% under any condition. 



With the depletion and desorption rates, defining fa to be the 
normalized fractional depletion, we can write an equation for fy, 

df d (r,t) 



X C0 (R, 1 = 0) 



10" ■ 



dt 



= A{T,t)f d {T,f)-^T,f){l-f d {T,t)), 



(6) 



in which we implicitly use the boundary condition that the sum 
of CO molecules in solid state and in the gas phase is constant. 

Equation [6] is a very general one and the rate functions can 
easily be ch a nged s o that the equation describes other reactions, 
van Weeren (2007) shows that in principle a whole set of such 
equations can be solved simultaneously u sing gas phase reac- 
tion rates from e.g., the UMIST database dWoodall et al]l2007l) 
although computational demand becomes an issue for a full net- 
work. 

3. Results 

To obtain abundance distributions we evaluate Eq. [6] using the 
gas temperature T g , the dust temperature T c /, and the mass den- 
sity (we assume that H2 is the main mass reservoir and that the 
gas-to-dust ratio is 100:1) from the hydrodynamical simulation. 
We assume that the cloud has been in hydrostatic equilibrium 
prior to collapse for 3 x 10 s yr (=1 x 10 13 s). The density profile 
before collapse is given by a spherical power-law with a uni- 
form temperature as described in Sect. [2] Starting with a pure 
gas phase abundance of 2 x 10~ 4 molecules relative to H2, we 
let the abundance evolve for 3 x 10 5 yr with a constant tempera- 
ture of 10 K and a static density profile. The new CO abundance 
serves as initial condition for the collapse, i.e., the start of the 
hydrodynamical simulation which defines t = 0. In the following 
we present three different ways to obtain the abundance profiles. 

3.1. Model 1: Drop abundance 

This approach is denoted the drop abundance model, because 
it follo ws directly from the method presented by J0rgensen et al.l 
(2005). In this approach we do not actually solve Eq.|6] but rather 
consider the rate equations [3] and |4] The rate equations are in- 
verted to obtain a characteristic time scale, and at any given time 
in the simulation, this time scale is evaluated based on the cur- 
rent density and temperature conditions. If the freeze-out time 
scale is shorter than the evaporation time scale, the gas phase 
abundance drops to a preset level (sufficiently low, i.e., several 
orders of magnitude) but if the actual age of the cloud is shorter 
than the freeze-out time scale, the gas phase abundance is kept 
at its initial value. Using this approach, the resulting profile is a 
step function where all CO is either in the gas phase or in the 
solid state. The movement of the trace particles does not affect 
the result since the abundance is determined only from the cur- 
rent local conditions. The effect of cosmic ray desorption is im- 
plicitly accounted for by not letting the abundance drop to zero, 
but only two or three orders of magnitude, although freeze-out 
conditions are met. 

The result is shown in Fig. [4] The top panel shows the ra- 
dial abundance profile along the disk mid-plane (z=0), while the 
lower panel shows the two dimensional distribution. The result 
of this drop abundance approach supports the idea of an evolu- 
tion of the CO abu ndances from a pre-ste llar core to a Class I 
object, proposed by Jeir gensen et al.l (|20 05). The maximum size 
of our depleted zone is ~6000 AU, just before the collapse sets 
in, and shrinks to ~500 AU (when averaged in radial and polar 
directions) as the cloud collapses and the disk is formed. The 
inner evaporated zone has a radius of ~1000 AU, in agreement 
with values derived in J0rgensen et al. (2005|). 
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Fig. 4. This figure shows the abundance distribution of Model 
1. The top panel shows the radial abundance profile through 
the mid-plane of the disk. The bottom panel shows the spatial 
abundance distribution. The shaded area is where molecules are 
frozen out. The unit on the y-axis of the top panel is fractional 
abundance. 



3.2. Model 2: Variable freeze-out 

The second method is similar to Model 1 except that instead of 
using the time scales, we use the actual rate equations which 
allows the trace particles to be in a mixed state. The values for 
A and £ are calculated for all r using the current n(r) and T(r) 
and Eq. [6] is solved to obtain the fraction of gas phase CO to 
solid state CO for each trace particle position. In this model, we 
still do not consider the history of the trace particles and thus the 
motion of the particles is irrelevant and the time dependence of 
£ and A in Eq. [6] drops out. Plots of the abundance distributions 
are shown in Fig. [3] The profiles of Model 2 are quite similar to 
the ones from Model 1, except for the gradual change from pure 
gas phase to pure solid state. Especially the inner edge is well 
approximated by the step function of Model 1 because of the 
exponential factor in the evaporation, described by Eq. @] The 
main difference between the two models is the outer edge, and 
especially at earlier times, where the gradient is more shallow. 




3.3. Model 3: Dynamical evolution 

Our final model takes the motion of the trace particles and their 
chemical history into account. Freeze-out and evaporation rates 
are calculated dynamically for each time step and composition 
of the particles is used as initial condition for Eq. [6] For exam- 
ple, the depletion in a certain area represented by a trace particle 
depends on where this particle has come from, which we know, 
because the trace particles follow the flow calculated in the hy- 
drodynamical simulation. Two trace particles may pass by the 
same region at the same time and experience the same desorp- 
tion and depletion conditions. However, they may have come 
from different areas, e.g., one may have spent time in a warmer 
region and the other may come from a colder region and there- 
fore they will leave the place with different compositions. 

Figure [6] shows the abundance distributions for this model. 
Surprisingly, the profiles are not very different from those of 
Model 2. A difference only occurs when the dynamical time 
scale is comparable to the chemical time scale. In this case, a gas 
phase trace particle may move out of a depletion zone faster than 
depletion can make any significant change to its composition, or 
vice versa. The evaporation front around the disk is also sharper 



compared to Model 2 which makes Model 3 appear almost like 
a step function in the latter part of the simulation. 

3.4. Cosmic ray desorption and increased binding energy 

We now investigate the effects of excluding the cosmic ray des- 
orption term (Eq. [5]) to see what effect this mechanism has on 
the abundances. We also increase the CO binding energy Et, in 
Eq. H]in order to emulate different ice mantle compositions. In 
the following we use Model 3 as a template. 

Cosmic ray desorption ensures that depletion can never reach 
100% as long as cosmic rays are able to penetrate the cloud and 
remove molecules from the ice matrix. If no other desorption 
mechanisms were present, dark cold molecular clouds would 
have a very high degree of depletion since therm al desorption 
is larg ely ineffective at low temperatures. However. ICaselli et al.l 
(1999) find no evidence for such high depletion factors and 
therefore some other desorption mechanism (such as for instance 
cosmic ray desorption) must play a certain role in these environ- 
ments. However, to see exactly how much evaporation cosmic 
rays account for, we have calculated Model 3 without the cos- 
mic ray desorption term. Figure [7] show these abundance dis- 
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Fig. 7. Abundance distribution of Model 3 without cosmic ray 
desorption and the default binding energy of 960 K. Otherwise 
similar to [4] 
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Fig. 8. Abundance distribution of Model 3, including cosmic 
rays, but with a CO binding energy of 1740 K. Otherwise similar 
tog] 



tributions. The profiles are largely similar to Model 3, but the 
amount of depletion is, as expected, about three orders of mag- 
nitude higher. 

Including cosmic ray desorption obviously has a big impact 
on the amount of gas phase molecules in the cloud, raising the 
gas-phase abundance significantly. However the effect of cos- 
mic ray desorption can be countered by a larger binding energy 
between molecule and grain surface. For example, raising the 
binding energy from 960 K to 1740 K almost entirely compen- 
sates the effect of cosmic ray desorption, when the depletion is 
averaged over the entire model. 

An increased binding energy applies when the CO ice man- 
tle is mixed with o ther species, typically water ice or NH3 
ice dFraser et alj |2004). Another way to obtain a spread in bind- 
ing energies is when the ice surfaces are not even, e.g., dif- 
ferently layered ices, varying mantle thickness, rough surfaces, 
etc. If the grains are more fractal, with cavities which can trap 
CO molecules, it will become harder to desorb these molecules, 
ef fectively raising the b inding energy. Laboratory experiments 
of ICollings etail d2004) even showed multiple peaks in the CO 
desorption spectrum, implying a range of binding energies. In 



Fig- Hit can be seen that the inner most evaporation zone is con- 
sistently and significantly smaller when the binding energy is 
increased compared to when cosmic rays are omitted, suggest- 
ing that the two effects can be disentangled by high-resolution 
observations. 

3.5. Comparison to observations 

We compare our results with observations of CO abundances and 
depletion factors. The fractional depletion fj of CO for the var- 
ious models are shown in Fig. [9] The highest depletion factors 
occur at the end of the pre-stellar core phase just before collapse. 
Observations show t hat depletion factors d ecrease with decreas- 
ing envelope mass ( J0r gensen et al.1 120021) . We indeed see the 
depletion decrease after the core begins to collapse due to CO 
evaporation close to the protostar. At the same time the envelope 
density decreases which also lowers the global fractional deple- 
tion, defined as the ratio of the total number of gas phase CO 
molecules to H2 molecules in the entire model. 

Freeze-out in the disk is most pronounced for Model 3 in- 
cluding cosmic ray desorption and a high CO binding energy. 
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Fig. 9. The global fractional depletion in the various models. The 
small bump in the curves around 0.5 t/f is due to a change in the 
luminosity output of the central source as the disk is formed. 



The same model, but with a binding energy of 960 K, gives a 
very low depletion factor at the start of t he collapse. Althoug h 
low values of fy have been measured (e.g., Bac mann et al.l2 002). 
most observations poin t toward a higher value, fj > 10, in 
pre-stellar cores (e.g., ICrapsi et al.l l2004t iBergin et all l2002t 
liessop & Ward-Thompson 1200 lb . The low value of /<? we find 
means that either the binding energy is larger than 960 K or 
that we ov erestimate the cosmic ray desorption rate. Indeed, 
IShen et all (120041) showed that the importance of cosmic ray 
heating has been overestimated by previous authors. A third pos- 
sibility is that CO can be converted into other species, resulting 
in a lower abundance, but since we only model CO we cannot 
quantify this effect. 

As mentioned above, a correlation has been found between 
the amount of depletion and envelope mass. The CO abundance 
drops when going from a pre-stellar core, to a Class 0, to a 
Cl ass I object. We comp are our results to the data presented 
by [jjirgensen et al. (2002) who derived global envelope abun- 
dances for 18 sources. Their result is shown in Fig.[l0] where 
open symbols are Class I objects, closed diamonds are Class 
objects, and the filled squares are pre-stellar cores. Superposed 
on the data are curves based on our models. To first order our 
models show the same trend. The observations show no evidence 
for freeze-out in the Class I objects but the number of this type 
of source is small and the derived envelope masses and abun- 
dances are somewhat uncertain. We also expect some amount of 
intrinsic scatter due to the unique environments and characteris- 
tics of each source. Again the drop abundance model and Model 
3 including cosmic rays and a low binding energy show too high 
CO abundance (or too little depletion). Unfortunately, our simu- 
lation begins with only one solar mass and therefore we do not 
catch the more massive Class and pre-stellar cores. 

We can also calculate gas column densities based on our 
abundance distribution models. Vertically integrated column 
densities, Nqo = fnco(z)dz, are shown at four different free- 
. [TT| The models differ mostly on the position 



fall times in Fig. 
of the evaporation front as a result of the different binding en- 
ergies used. The column de nsity profiles at t = can be com - 
pared to those presented by lAikawa et all (1200 ll I2003L l2005h . 
Although our underlying density model is different, the column 
densitie s are compar able, differing by an order of magnitude 
or less. Aika waet al.l show flat profiles toward the core center 



Fig. 10. In this plot we show the abundance as function of en- 
velo pe mass. Also shown in thi s plot are data points taken 
from lj0rgensen et al.l (120021 12005). The open symbols are Class 
I objects, the filled diamonds are Class sources, and the filled 
squares are pre-stellar cores. 
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Fig. 11. CO gas column densities of the models discussed in this 
paper. These profiles are seen through a pencil beam, and need 
to be convolved with an appropriate beam in order to be directly 
compared to observations. 



or slight drops within ~5000 AU, which is the same trend we 
find with our models. The only exception is our drop abundance 
model (Model 1) which gives an increasing column density to- 
ward the center. The drop occurs at -6000 AU, and within a 
6000 AU radius we have a constant abundance. The «h 2 density 
decreases quadratically with radius which means that the col- 
umn density of CO increases toward the core center. Th is is in 
disagreement with observations (e.g., Tafalla et al. 2004) where 
flat or decreasing integrated intensity profiles are found in core 
centers, which can be interpreted as a drop in the column density. 

Column densities drop around three orders of magnitude 
from th e center toward the f reeze-out zone in the disk, compa- 
rable to lAikawa et al.l ([1996), although column densities in our 
results are generally more than one order of magnitude higher. 
This is simply explained by the differences in the adopted disk 
masses. A typical value adopted for disk masses is 10~ 2 M , 
while the disk in our simulations grows to an unrealistically high 
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Fig. 12. Same as Fig.QT|but for CO ice column densities. 

mass of about 0.4 M Q . Also our column density distributions 
are somewhat more complicated as we do not use a s i mple ana- 
lytic description of a non-evolving disk. lAikawa et al.l (1 19961) get 
drops in column density due to freeze-out in the disk around 200 
- 300 AU, similar to our results of 100 - 500 AU, depending on 
which model we use. 

In a similar way we can calculate the column densities for 
the CO bound to ice. Plots of the ice column densities are shown 
in Fig. [12] Before collapse the column density of CO ice increase 
inwards due to the shorter freeze-out timescales there. When the 
central temperature rises, the column density flattens and drops, 
while in the protoplanetary disk the column density sharply in- 
creases. Overall, the column density drops in the envelope be- 
cause of the decreasing density (which also results in a lower 
freeze-out rate). It is difficult to compare these results with obser- 
vations since often cold (unrelated) clouds, where CO is frozen 
out, lie in front of protostars. Considering this, observed column 
densities of CO ice lie often around 10 18 cirT 2 in reasonable 
agreement with our results. 

3.6. Emission lines 

Simulated spectral lines are a direct way to test our models 
since these can be compared directly to observed spectra. As 
mentioned in Section [2] we use the hydrodynamical time snap- 
shots including the abundance distribution mapped onto the orig- 
inal grid as input models for our 2D molecular excitation code 
RATRAN. This code solves the level population using an accel- 
erated Monte Carlo method, which is then ray-traced to make 
frequency dependent intensity maps. After beam convolution, 
spectral profiles can be extracted from these and directly com- 
pared to observations. We add a mean field of 0.2 km s _1 to the 
snapshots during the radiation transfer calculations to emulate 
the presence of micro-turbulence. 

An in-depth analysis of spectra based on the hydrodynami- 
cal simulation is given in a separate paper (Brinch et al. in prep.) 
and only example spectra of the abundance models derived in 
this paper are presented here. Figure [T3] shows the CO 7=3- 
2 transition of the four time snap shots which have been used 
in the previous figures. The top three panels correspond to the 
three abundance models including cosmic ray desorption, while 
the lower three panels show spectra based on model 3 with- 
out cosmic ray desorption and with increasing binding energies. 
Figure [14] has a similar layout, but for the less optically thick 
isotopologue C ls O. For comparison, CO and C ls O spectra, cal- 
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Fig. 13. Time series of CO 3-2 spectra. The top three panels 
show the three different abundance models. The three lower pan- 
els show model 3, including cosmic ray desorption using three 
different binding energies. Time is shown by increasing line 
thickness and the four lines correspond to the four adopted time 
snapshots (0.0, 0.5, 1.5, and 2.5 iff). Each line has been offset 
by 0.5 K. 
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culated with a constant abundance profile are shown in Fig. [15] 
All spectra are seen under a 90° inclination (edge-on; similar to 
the orientation of the model in Fig. [TJ, and they have all been 
convolved with a 10" beam which is a typical value for single- 
dish sub-millimeter telescopes. A distance of 140 pc (Taurus star 
formi ng region) is adopted . We use an isotopic 16 0: I8 ratio of 
560 (I Wilson & Roodll 1994 to calculate all C ls O spectra. 

First of all, we see that the spectra in the top row of Fig. [13] 
are largely similar to each other and also to the constant abun- 
dance profiles in Fig. [15] The reason is that CO gets optically 
thick before the abundance goes down and thus they are inde- 
pendent of the underlying abundance model. Therefore we can- 
not use CO spectra to distinguish between the three models and 
we cannot determine the depletion factor using CO. Excluding 
cosmic rays or increasing the binding energy have the effect on 
the spectra of lowering the intensity, and in this case CO can be 
used to constrain the models. 

The line profile of the earliest snapshot (the thinnest line) 
is seen to have a very boxy shape with a flat line center. This 
is entirely due to the hydrodynamical model, which has not yet 
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Fig. 15. Time series of CO and C ls O spectra with a constant 
abundance profile. Time is represented by increasing line thick- 
ness. 



evolved from its initial isothermal state in the region where the 
line get optically thick. Basically we probe the same temperature 
through a range of velocities, which gives rise to the fiat top. The 
line of the earliest snapshot is also considerably more narrow 
than the other lines. This is because the velocities in the model 
have not yet had the time to build up. 

The C ls O spectra however are seen to differ considerably 
from the ones calculated from a constant abundance model. They 
also differ, although somewhat less, from each other when the 
three abundance models are compared. The reason is that the 
lines are optically thin and therefore a drop in the abundance has 
a direct impact on the spectral line, no matter where the drop 
occurs. 

The spectra shown here are all simulating single-dish mea- 
surement, where the emission of the source has been smeared 
out in a relatively large beam. If instead an interferometer is 
used, beam sizes better than 1" can be obtained, which allows 
us to map the line-of-sight column densities. With such observa- 
tions, it should in principle be possible to map out the abundance 
distributions of optically thin species and constrain the models. 

4. Discussion 

The abundance profiles presented here depend on the adopted 
hydrodynamical model. While we believe that the overall be- 
havior of our results are generally applicable, details such as the 
exact timescale and the absolute values of the fractional deple- 
tion depend on the initial conditions of our simulation such as 
mass and angular momentum and also on the particular hydro- 
dynamical scheme which we have used. 

A main concern regarding the hydrodynamics is the inter- 
nal radiation transfer module that is used to calculate the tem- 
perature structure of the simulation. This is done with an ap- 
proximate method which is known to overestimate the dust 
temperature by a few Kelvins compared to calculation done 
by a more accurate continuum radiation transfer code (Visser, 
priv. comm.). While a few Kelvin have little or no effect on the 
hydrodynamics and therefore on the evolution of the cloud, it can 
certainly make a difference for the evaporation front which has 
a strong temperature dependence. While this will affect the de- 
tailed behavior of the abundances, it will not change the general 
trends discussed in this paper. 



Although in this paper we only concerned ourselves with the 
CO gas-phase abundance, as mentioned in Sect. [2] our method 
can be easily extended to include other species and even surface 
reactions. However, running the full chemical network is con- 
siderably more computationally demanding. Whereas the CO 
model runs in a few hours on a standard desktop computer using 
9 x 10 5 trace particles, a test run of the full chemical network 
showed that it takes a factor of 100 longer to run using a factor 
of 1000 less trace particles. Nevertheless, such large scale chem- 
ical models present a powerful way to obtain realistic abundance 
profiles for hydrodynamical simulations which instruments like 
the Atacama Large Millimeter Array (ALMA) can test at high 
spatial resolution. 



5. Conclusion 

In this paper we have investigated different scenarios for the de- 
pletion of CO in a collapsing cloud. We have used a hydrody- 
namical code to describe the dynamical evolution of the cloud 
and the chemistry has been solved using a trace particle ap- 
proach. 

Of the various models shown in Sect. [3] the one which ex- 
clude cosmic ray desorption and the ones using higher bind- 
ing energies differ the most. The drop abundance model pro- 
vides a good approximation to Model 2 and 3. Therefore, if 
it is possible to fit a drop profile to observations, as was done 
by |j0rgensen et al.1 d2005l) . it should be possible to run a hydro- 
dynamical simulation with customized initial conditions, so that 
it reproduces this profile and thereby putting an observational 
constraint on the hydrodynamics. Also in cases where an analyt- 
ical model is used, using the drop abundance approximation is 
quite good. 

We furthermore conclude that the freeze-out chemistry has 
little effect on optically thick spectral lines, especially when ob- 
served in low resolution. It is therefore not entirely unreasonable 
to model such lines using a constant abundance model disregard- 
ing the chemistry. However, one should be careful when using 
optically thin lines to determine which average gas phase abun- 
dance to use for a constant abundance model, since the average 
abundances derived from such lines are indeed depending on the 
freeze-out chemistry. The maximum discrepancy between the 
flat abundance profile model (Fig.[T5ll and the freeze-out models 
shown in Fig. [14] is a factor of three in the integrated intensity. 
Such a change in the average abundance does not affect the opti- 
cally thick lines a lot, so depending on the accuracy needed, us- 
ing a flat abundance distribution to model data in which freeze- 
out is present may not be entirely unreasonable. 

Finally, the effect of the chemical depletion is considerably 
enhanced when the resolution of the observations is increased, 
making it more important to get the abundance profiles right. 
This will especially be important for future ALMA observations 
(with a resolution of ~ 0.01") where the level of detail will be 
so great that getting the exact abundances will be crucial for 
interpreting the data. 
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